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Abstract 


A novel approach is proposed for high-fidelity modeling of progressive damage and failure in 
composite materials that combines the Floating Node Method (FNM) and the Virtual Crack 
Closure Technique (VCCT) to represent multiple interacting failure mechanisms in a mesh- 
independent fashion. In this study, the approach is applied to the modeling of delamination 
migration in cross-ply tape laminates. Delamination, matrix cracking, and migration are all 
modeled using fracture mechanics based failure and migration criteria. The methodology 
proposed shows very good qualitative and quantitative agreement with experiments. 

1 Introduction 

Damage in composite materials generally occurs as a combination of different and interacting 
failure mechanisms, e.g. delamination and matrix cracking. Capturing these interactions 
accurately is essential to confidently model and predict progressive damage and failure. Many 
approaches have recently been proposed that explicitly model different failure mechanisms and 
attempt to capture their interaction [1-3]. However, no single strategy has yet surfaced as totally 
satisfactory. The present approach combines the Floating Node Method (FNM) [4] with the 
Virtual Crack Closure Technique (VCCT) [5,6] to explicitly account for different failure 
mechanisms and their interaction. Delamination, matrix cracking, and migration events are all 
modeled with the same FNM element using fracture-mechanics-based failure and migration 
criteria. The approach is validated using recent experimental results in which a setup capable of 
isolating a single complete migration event was developed [7]. Migration is defined as the 
transition of a delamination from its current ply interface to a new interface via transverse ply 
cracking. The results obtained offer a simple, yet very challenging, benchmark for the validation 
of numerical and analytical approaches aimed at high-fidelity modeling of composite materials, 
such as the one proposed in the present paper. 

2 Floating Node Method (FNM) 

2.1 Introduction 

The FNM is a recently proposed numerical method, capable of representing multiple evolving 
discontinuities in solids [4], Compared to other methods used for the same purpose, such as 
extended Finite Element Method (XFEM) or the Phantom Node Method (PNM), it provides the 
following key advantages. First, FNM does not introduce an error on the crack geometry when 
mapping to natural coordinates. Second, it does not require numerical integration over only part 
of a domain. Third, it is ideally suited for the representation and incorporation of multiple and 
complex networks of weak, strong and cohesive, discontinuities. Fourth, FNM yields the same 
solution as a finite element mesh where the discontinuity is represented explicitly, and finally it 
is conceptually simpler than the PNM or XFEM. 
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2.2 Element formulation 


The static equilibrium of a body with volume ft under body forces with density f (acting on ft) 
and traction t acting on the boundary can be expressed in the weak form as: 

J n e T (v)a(u)dft = J n v T fdft + v T tdr n (1) 

where u is the displacement vector; v is the test function; e is the strain tensor related to u 
through the differential operator relative to Cartesian coordinates £ x as € = £ x (u); and a is the 
stress tensor related to the strains through Hooke’s law as a = De, with D being the constitutive 
tensor. In the Floating Node Method, each real node “i” is characterized by its nodal coordinates 
Xj and associated Degrees of Freedom (DoF) qj. In addition, an FNM element contains a suitable 
number of floating DoF without pre-defined associated nodal position vectors. These additional 
floating nodes are used to represent discontinuities. Their number varies with the number and 
type of discontinuities, weak or strong, modeled within each element. A strong discontinuity is 
defined as a jump in a field quantity (e.g. displacements), while a weak discontinuity is defined 
as a jump in the gradient of a field quantity (e.g. strains). Figure 1 shows an example of such an 
element with four nodes and four additional sets of floating DoF required to represent a strong 
discontinuity. 

2.2.1 Element formulation without weak/strong discontinuity 


If there is no discontinuity to be modeled by the element (e.g. before failure), the formulation is 
the same as in the standard Finite Element Method (FEM). The vector of nodal coordinates is 
given as x n . In the case of Figure 1, x is given by: 


Xn = [X!,X 2 ,X3,X 4 ] 

The Jacobian of the transformation from physical coordinates x to natural coordinates ^ is: 
dx dN 


J = dTd^ 


( 2 ) 

( 3 ) 


where N is a matrix of shape functions. Assuming an isoparametric formulation, the 
displacement field u is related to the real DoF q through: 

u = Nq (4) 


2 


In the case of Figure 1, q r = [q 1( q 2 , q 3 , q 4 ]. The strains can also be expressed in terms of the 
real DoF q as: 

€ = £ x (u) = X^(N)J- 1 q = Bq (5) 

where B = £^(N)J -1 , leading to the stiffness matrix: 

K = J B t D Bdet(J)dS (6) 

where S is the integration domain (in natural coordinates) corresponding to ft (in physical 
coordinates), and to the vector of nodal forces: 

Q = J N T f det(J)dS + J N T fdet(J)dS (7) 

where T s is similarly the boundary corresponding to The weak form of equilibrium, Equation 
1, becomes: 


Kq = Q 


( 8 ) 


2.2.2 Element formulation with weak/strong discontinuity 

Once a discontinuity in the element is defined, the element is split in two or more partitions 
(depending on the discontinuity). Without loss of generality, a case in which the element is split 
in two partitions, ft A and ft B , is illustrated in Figure 1. For each partition, ft A and ft B , a vector of 
nodal coordinates, x nA and x nB , is defined. For the case in Figure 1: 


x n A = [x5,x 6 ,x 3 ,x 4 ] 

(9) 

x n B = [x 1( X 2 ,x 6 ,x 5 ] 

(10) 

Each partition has a separate Jacobian: 


dx dN 


,A = d^ = d^ XnA 

(11) 

dx dN 


,B = d^ = d^ XnB 

(12) 


The displacements u A and u B , in partitions ft A and ft B , respectively, are interpolated separately 
from the DoFs q A and q B associated with each partition: 


u A = Nq A and u A = Nq A 


(13) 


These DoFs, q A and q B , did not have an associated position beforehand, and therefore were 
considered to be “floating”. As the discontinuity is defined, they are linked to the respective 
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position vector. In the case of Figure 1, which represents a strong discontinuity, 
q A = [cj5, cj6> *14] and q B = [qi, Q2- q6/« Qs/l- Note that there are four different sets of 
floating DoF: q 5 is different from q 5 , and q 6 is different from q 6 ,. If a weak discontinuity were 
to be modeled, only two sets of floating DoF would be included in the element, and the DoF with 
a prime would coincide with those without a prime. 

The strains then become: 


€ a — Ac( u a) — £^(N)Ja 1( 1a — B A q A 

(14) 

€ b = £x( u b) = ^(N)J b x q B = B B q B 

(15) 

The stiffness matrices for partitions fi A and ft B are: 


K a = BlDB A det(J A )dS 

(16) 

K b = B B DB B det(J B )dS 

(17) 

and the force vectors: 


Qa = J s N T f det(J A )dS + / rs N T f det(J A )dS 

(18) 

Qb = f E N T f det(J B )dS + J rw N T f det(J B )dS 

(19) 


Finally, the equations of equilibrium can be written as: 


KaQa — Qa ar *d K B q B — Q B 


( 20 ) 


2.2.3 Element topology and assembly 

To illustrate this approach, a floating node element capable of modeling two weak interfaces 
with an arbitrary crack between them was implemented, as shown in Figure 2(a). The nodal 
position of the initial integration domain, and the real and floating DoF are indicated in the 
figure. Each floating DoF is associated with either an edge or the inner domain of the element. 
Once weak/strong discontinuities are detected, the floating DoF at the edges are used to 
determine the solution for the nodal positions given by the intersection of these discontinuities 
with the element edges (Figs. 2(b) to 2(e)). The additional inner floating DoF are used to 
determine the solution for the nodal positions originated by the intersection of multiple cracks 
within the element (Fig. 2(f)). In Figures 2(b) to (f) only the floating DoF used are represented. 
In each case, the floating DoF that are not used have no assigned position, only a topological 
relation to an edge or to the inner domain. During assembly, inner floating DoF can be removed 
from the analysis through static condensation. The floating DoF associated with the edges are 
assembled with the corresponding DoF of neighboring elements, and will therefore have a 


4 


unique position in the global DoF vector. All floating DoF are then used as required to model 
weak/strong discontinuities as they evolve throughout the simulation. 

The approach can be extended to model an arbitrary number of weak and strong discontinuities 
within an element, and therefore an arbitrary number of plies and interacting cracks, provided the 
topology is adequately defined, and a sufficient number of floating DoF are used. 

3 Virtual Crack Closure Technique (VCCT) 

In the present work, the FNM method has been coupled to the Virtual Crack Closure Technique 
(VCCT) [5]. An extensive review of VCCT is provided in [6]. Traditionally, VCCT is used to 
obtain energy release rates in cases where the crack path is known beforehand and can be aligned 
with the elements’ edges, such as the case for de lamination. 

In VCCT, the Mode I and II energy release rates are obtained from [6]: 



where F n and F t are the normal and tangential components of internal force vector F at the crack 
tip; IQvI and are the normal and tangential components of the displacement jump [q] 
between the nodes immediately behind the crack tip; and a x and a 2 are the lengths of the crack 
in the elements behind and ahead of the crack tip, respectively (Fig. 3). When combining VCCT 
with FNM, the forces and displacements needed to obtain the energy release rates, Eqs (21) and 
(22), are computed at the floating DoF as the crack develops (Fig. 3). Additionally, floating DoF 
are also added along a virtual crack plane up to a distance r, named the ‘enrichment radius’, r enr . 
The associated floating nodes provide additional degrees of freedom along this path, which 
enable the accurate modeling of the deformations near the crack tip [4], At the end of the 
enrichment radius, the displacements of the floating DoF are interpolated from the displacements 
of the real DoF. In the present work, the enrichment radius was chosen to be equal to the largest 
in-plane dimensions of the numerical model. 

4. Delamination propagation 

In this paper, delamination is modeled using VCCT to compute the energy release rates G, and 
G n at each delamination front position. These are then used in a failure criterion: 

Gn) = — 1 = 0 (23) 

' J c 

where G T = G, + G n , and G c is the critical energy release rate given by [8]: 

G c = G lc + ( G Uc — G lc ) ^ (24) 
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For delamination growth according to this criterion, G c is assumed to equal the critical energy 
release rate of the interface, G l c nt . The delamination front is advanced by one element along the 
projected crack path when G T = Gc nt . 


5 Matrix crack propagation 

In the present study, matrix crack propagation is also modeled with YCCT. However, unlike 
delamination, the matrix crack is assumed to propagate following a Mode I path in the through- 
thickness direction, as supported by experimental evidence [7]. Therefore, the failure criterion is 
assumed to be given by: 

/(G,) = y-- 1 = 0 (25) 

where the critical energy release rate is given by G c = G, c . 

At each crack growth increment, the Mode I crack path orientation is approximately determined 
using a maximum tangential stress criterion. This criterion, typically written in terms of the 
Mode I and II stress intensity factors [9], can also be written in terms of energy release rates G, 
and G n , where the angle 0 that maximizes the tangential stress Ogg is obtained by evaluating: 



and for the two solutions obtained choosing the one that maximizes Ogg given by: 

000 = yfGi {3 cos (0 + cos (y)} + yG^Sgn(F t ) {-3 sin (0 - 3 sin (y)} ( 27 ) 


where F t is the tangential component of the internal force vector F at the crack tip. The tangential 
stress a ee relates to Og'g by: 


_ 1 E* _ 
where, 


(29) 


E* 



plane stress 
plane strain 


(30) 


E is the Young’s modulus; and v the Poisson’s ratio. Once the angle is determined, the criterion 
given by Eq. (25) is assessed. If the criterion is met, the crack advances by one element along the 
projected crack path. 
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6 Migration 

6.1 Delamination to matrix crack - migration criterion 


Extensive research has been done on migration of cracks from a bi-material interface. He and 
Hutchinson [10] have proposed a necessary condition for migration to occur: 


r B 

u max ^ U T 

~Gf 


Gn nt 


G' c nt 


(31) 


This condition states that for a crack to migrate from an interface into a given Material B, as 
shown in Figure 4, the ratio between the energy release rate associated with its growth into that 
material G^ ax , (at an angle oj max that maximizes its value) and the critical energy release rate of 
that material G B , has to be higher than the ratio between the energy release rate associated with 
the growth of the crack along the interface G| nt , and the critical energy release rate G ! c nt of that 
interface. 

With the approach proposed, evaluating Eq. (31) would require more than one FE analysis at 
each interface crack position. In each of these analyses, a branched crack would need to be 
initiated at different tentative angles co in order to determine the angle, oj max , that would 
maximize G B . 

In the present work, an alternative criterion is proposed. With this criterion, all variables needed 
to determine whether migration occurs, can be directly obtained from the FNM results at the 
each interface crack position. The criterion states that: (i) the sign of the tangential component of 
the internal force vector F t , associated with Mode II shear displacement, dictates the direction 
and therefore the material into which the interface crack migrates and (ii) migration can only be 
completed if it is energetically favorable. The criterion is given by Mg : {Gf nt , F t } -> {—1,0,1} as: 


M§ = 



Sgn(F t ), 


Sgn(F t ), 


Sgn (F t ) < 0 


Sgn (F t ) > 0 


(32) 


If Mg = 0 the interface crack does not migrate, and if Mg = — 1 or Mg = 1 the crack migrates 
into material A or B, respectively. In the present form, Eq. (32) assumes that internal force vector 
is computed at the lower surface (Material B). If the internal force vector is computed at the 
upper surface the signs of the inequalities in Eq. (32) need to be reversed. Depending on the sign 
of F t , the ratio between the energy release rate determined at the interface, and the critical energy 

Qlnt Qlnt 

release rate of Material A or B, or is compared to the ratio between the energy release 

G c G c 

Qlnt 

rate determined at the interface and the critical energy release rate of the interface, If the 

g c 


1 


Qlnt Qlnt Qlnt 

ratio or is greater than -f^, the necessary condition for migration is met, and Mg = — 1 

or Mg = 1. With this criterion, all variables needed to determine whether migration occurs, F t 
and Gj lt , can be directly obtained from the FNM results. Note that although this criterion is 
necessary, it is not sufficient for migration to occur unless the ratio being assessed is greater than 
one. 

In the present work, the interfacial crack, or delamination, propagates between a 0° and a 90° 
ply. Assuming Material A is the 0° ply, and Material B is the 90° ply, the critical energy release 
rate to propagate a crack into a 0° ply is much greater than the critical energy release rate of the 
interface, since fiber fracture is required for that process to occur, i.e. G°° = G? » G l c nt [11]. 

Qlnt Qlnt 

Thus, is always smaller than -j^. Therefore, even if the sign of F t favors migration into 

Material A, migration will not occur, since it is not energetically favorable. Typically for 
composite laminates, the critical energy release rate of an interface crack (delamination) and of a 
matrix crack are of the same order of magnitude, i.e. Gc°° ~ Gc nt [12]. Hence, Gc nt is assumed 
to be given by Eq. (24) and Gc°° to be equal to G Ic . Therefore, assuming Material B is the 90° 
ply, provided the sign of F t favors migration into Material B and it is energetically favorable, 
migration can occur. 

6.2 Matrix crack to delamination - migration criterion 

When a matrix crack reaches a weak interface, it is assumed to trigger delamination. This 
delamination is at first contained within the FNM element, as illustrated in Figure 2(f). In the 
following steps, the delamination will propagate, or migrate, as detailed in Sections 4 and 6.1, 
respectively. 

7 Simulation of delamination migration specimens 

7.1 Delamination migration test configuration 

In reference [7] a test aimed at investigating delamination migration was proposed. Its 
configuration is illustrated in Figure 5. The specimen possesses three main features that permit 
the controlled observation of delamination growth followed by migration to another ply 
interface. First, the specimen geometry is in the form of a beam with the intent of promoting 
uniform delamination growth and migration across the specimen width. Second, the specimen 
contains a PTFE insert (acting as an artificial delamination) at an interface between a 0° ply 
(specimen span direction) and a stack of four 90° plies (specimen width direction) (Fig. 5). This 
provides an opportunity for the delamination to migrate to another ply interface by kinking 
through the 90° ply stack. Third, the specimen can be loaded in a manner to cause delamination 
growth from the PTFE insert that eventually migrates to another ply interface. This sequence of 
fracture events is made possible by the way in which specimen loading affects shear stresses 
acting across the delamination front. 
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7.2 Finite Element Model 


The experiments presented in [7] were simulated using the finite element solver 
ABAQUS/Standard® (Implicit). Plane-stress and a small-displacement formulation were used. 
The floating node element (Fig. 2) was developed and implemented as a user defined element 
(UEL) in ABAQUS, and applied in the center region of the model, as shown in Figure 6. The 
specimen layup used in the experiments and in the model, is also given in Figure 6. In the layup 
description T represents the PTFE insert. Each block of plies of the same orientation was 
modeled with a separate element CPS4 (through-thickness), except for the highlighted plies at 
the center of the model (Fig. 6). These were modeled within a single FNM element, as illustrated 
in Figure 2. The material properties used are given in Tables 1 and 2. A prescribed displacement 
was applied at different positions, L, at a single node, to simulate the displacement-controlled 
test. Four load positions were simulated corresponding to L = (a 0 , l.la 0 , 1.2a 0 , 1.3a 0 ] (Fig. 6). 
In the experiments, the specimens were subject to moisture intake from the ambient air leading to 
relaxation of residual thermal stress. Therefore, residual thermal stresses were not considered. 

7.3 Load-displacement response and sequence of events 

The computed load-displacement curves obtained for different values of L are given in Figure 7. 
Even accounting for the compliance of the test fixture, measured in [7], all load-displacement 
curves obtained from the simulations show a slightly stiffer response than what was observed 
experimentally. The maximum load is in general very well predicted for the cases L = a 0 to 
L = 1.2a 0 , and slightly over-predicted for the case L = 1.3a 0 (Figs. 7(a) to 7(d)). For L = a 0 , 
the simulations predict unstable delamination after the peak load (Fig. 7(a)). This unstable 
delamination stops just before migration. After the migration event, delamination continues 
stably along the next 0°/90° interface. A similar sequence of events was observed 
experimentally [7]. However, in the experiments, after the first unstable event, further loading 
was necessary before migration occurred. For the cases L = 1.1 a 0 to L = 1.3a 0 , the simulations 
predict that (i) the first peak in load is followed by a region of stable delamination growth (Figs 
7(b) to 7(d)) that increases with the load offset, and that (ii) stable de lamination is followed by a 
region of unstable delamination, which corresponds to the load-drop in Figures 7(b) to 7(d). The 
migration event occurs within this load-drop. After the migration, unstable delamination 
continues propagating along the next 0°/90° interface. Loading the specimen further eventually 
leads to a transition from unstable to stable delamination propagation. This sequence of events is 
very similar to what was reported experimentally [7]. The main difference is that in the 
experiments, for the cases L = 1.2a 0 and L = 1.3a 0 , after the peak load, delamination 
propagates through a series of unstable events, rather than stably, as predicted by the simulations. 

7.4 Migration location 

Figure 8(a) provides an illustration of the typical crack paths predicted by the simulations and 
obtained experimentally. Figure 8(b) compares numerical and experimental results for the 
distance from the load-application point at which migration onset occurs, A K , as a function of the 
normalized load offset, L/a 0 . Overall, except for L = a 0 , migration is predicted to occur closer 
to the load application point than what was observed experimentally. The simulations also 
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predict a decreasing trend in A K with L/a 0 , which is similar to what was observed 
experimentally for the cases L = l.la 0 to L = 1.3a 0 . 

7.5 Matrix crack path 

In Figure 9, the simulated matrix crack paths obtained for the different load positions was 
superimposed on an image of a typical crack path observed experimentally, where L = a 0 . The 
predicted matrix crack path agrees well the crack path observed experimentally [7]. The 
simulations do not predict variation in the matrix crack path with a change in load offset. 

7.6. Discussion 

All load-displacement curves obtained numerically show a slightly stiffer response than the 
experimental results. In the simulations, the engineering properties used were obtained from 
uniaxial test data. These are known to be higher than their flexural counterparts [13]. Since the 
specimen is predominantly loaded in flexure, the difference between uniaxial and flexural 
properties is therefore likely to be the cause for the overall stiffer response obtained in the 
simulations. Overall, the simulated load-displacement curves, sequence of events, and failure 
morphology are in very good agreement with the experimental results. For the case L = a 0 , 
experimental results show an increase in load at the end of the first unstable event, before 
migration. This increase in load is likely to be caused by the development of fiber bridging 
during the unstable part of the delamination, which was not accounted for in the model. For the 
load cases L = 1.2a 0 to L = 1.3a 0 , the only significant difference between simulations and 
experiments is that delamination growth onset and propagation, prior to the final load-drop see 
(Figs. 7(c) and 7(d)), does not occur stably as predicted in the simulations, but rather through a 
series of unstable events. In these experiments the unstable delamination growth onset was likely 
due to the presence of resin rich pockets at the tip of the PTFE insert, while the subsequent 
unstable events, prior to the final load-drop, may be due to features accompanying delamination 
growth, such as fiber bridging. Both of these were not accounted for in the simulations. 

The migration location was slightly under-predicted for all cases except L = a 0 . Nevertheless, 
the migration location trend observed experimentally, for the cases L = 1.1 a 0 to L = 1.3a 0 , was 
well captured. This further supports the validity of the approach, and in particular, of the 
migration criterion proposed. The latter is additionally confirmed by the good agreement 
between the predicted matrix crack path and that observed experimentally. This also suggests 
that the criterion used to determine the matrix crack orientation (as it propagates in the through- 
thickness direction) is adequate. 

8. Conclusions 

A novel approach that combines the FNM (Floating Node Method) with VCCT (Virtual Crack 
Closure Technique) has been demonstrated to adequately capture delamination migration in 
cross-ply laminates. Delamination, matrix cracking, and migration, are all modeled with the 
same FNM element using fracture mechanics based failure and migration criteria. Results show 
very good agreement with the maximum load, and migration location obtained experimentally. 
In addition, simulations capture the overall failure morphology and typical crack path, observed 
in the experiments. 
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Table 1. Elastic properties, IM7-8552 [14]. 


E n 
(GPa) 

^22 ~ E 33 
(GPa) 

V 12 ~ V 13 V 23 Gl2 ~ Gl3 

(GPa) 

£*23 

(GPa) 

a ± 

(K- 1 ) 

a 2 = a 3 
(K- 1 ) 

161.0 

11.38 

0.32 0.44 5.17 

3.98 

lxlO -5 

31xl0 -6 


Table 2, Critical energy release rates for delamination of IM7-8552 [15]. 


Gic 

Gnc V 

(kj/m 2 ) 

(kj/m 2 ) 

0.21 

0.77 2.1 


Without discontinuity With weak/strong discontinuity 




£)A 


Figure 1 . Overview of the Floating Node Method (FNM). 
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(a) topology of the floating node element (b) two weak interfaces, e.g. due (c) delamination 

to plies with different 


orientations 



(d) delamination with migration to matrix crack 


(e) matrix crack 


(f) matrix crack with 
migration to delamination 


Figure 2. Floating node element used showing different weak/strong discontinuities scenarios 



Figure 3. Virtual Crack Closure Technique and Floating Node Method for arbitrary crack propagation. 



Figure 4. Bi-material interface. 
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Applied 

load 


All units in mm. Not drawn to scale 
Specimen width, B=12.7mm 
90°- specimen width direction 
0°-speciment span direction 



////////////// 


Figure 5. Schematic of delamination migration test configuration [7]. 
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5.5 mm 


Stacking 

[90 4 /0 3 /(90/0) 2s /0 2 /0/90 4 / T /0/90 4 /0/0/(90/0) 2s /0/0/90 3 /0/90] | 

Plies modeled within the 
FNM element 




90° 


^ 0 


Figure 6. Finite element model, and boundary conditions of migration test specimen [7]. 
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(a) load offset: L = a 0 , (b) load offset: L = l.la 0 , 



(c) load offset: L = 1.2 a 0 , (d) load offset: L = 1.3a 0 , 

Figure 7. Experimental and predicted load-displacement plots. 
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♦Experiments 


1.1 1.2 1.3 


Load offset, LI a 0 

(a) qualitative comparison between simulated and (b) variation of A K with normalized load offset L/ a 0 

experimental crack paths 

Figure 8. Crack path and migration location. 



Figure 9. Comparison of experimentally observed crack path for L = a 0 , and simulated crack 
paths for different load-offsets, L. 
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